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Abstract 


From the overlap lattice quark propagator calculated in the Landau 
gauge, we determine the quark chiral condensate by fitting operator product 
expansion formulas to the lattice data. The quark propagators are computed 
on domain wall fermion configurations generated by the RBC-UKQCD Col- 
laborations with Ny; = 2 + 1 flavors. Three ensembles with different light 
sea quark masses are used at one lattice spacing 1/a = 1.75(4) GeV. We 
obtain (uw)M9(2 GeV) = (—305(15)(21) MeV)? in the SU(2) chiral limit. 
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1 Introduction 


The strong interactions among quarks and gluons have two prominent features at 
low energies: confinement and chiral symmetry breaking. The quark chiral con- 
densate (yw), which is in the light quark massless limit, is the order parameter of 
the spontaneous chiral symmetry breaking in Quantum Chromodynamics (QCD), 
the theory describing strong interaction. Furthermore X = —(ww) is one of the 
two low energy constants of chiral perturbation theory, the low energy effective 
theory of QCD, at leading order. The quark chiral condensate also appears in 
QCD sum rules and is an important input parameter. 

Thus there have been many determinations of the chiral condensate from dif- 
ferent ways by using lattice QCD, which is the nonperturbative method to solve 
QCD from first principles. See, for examples, Refs (4] B (6) [71 [S] [9]. A 
review of the evaluations of the chiral condensate on the lattice can be found in 
Ref. [10]. 

In this work, we determine the SU(2) low energy constant X by comparing 
the Operator Product Expansion (OPE) of the quark propagator in momentum 
space in the continuum MS scheme with the lattice calculation of the propagator in 
Landau gauge. This strategy was used by the ETM Collaboration in a calculation 
with two flavors of dynamical Wilson twisted mass fermions [11]. Our analysis 
is based on 24-1-flavor domain wall fermion configurations and overlap valence 
quarks. There were also analysis using the staggered fermions [12], the OPE of 
the pseudoscalar vertex [14] and the OPE of the quark propagator in coordinate 
space [15]. 

Our final result obtained at one lattice spacing is 3!/? = 305(15)(21) MeV in 
the MS scheme at the renormalization scale 2 GeV. Here the first error contains 
uncertainties from statistics, the lattice spacing and truncation effects in perturba- 
tive calculations. The second error is an estimation of the O(a?g?) lattice artifacts 
in our data. 

In the rest of the paper, we first discuss the OPE of the quark propagator in 
the MS scheme in Sec. 2| Then our lattice setup is given in Sec. |3| The analysis 
of the quark propagator and the results of the chiral condensate are presented in 
Sec. |4| Finally we summarize in Sec. 


2 OPE of the quark propagator 


For a quark field with mass m,, its propagator in momentum space S,(p) can 
be written as 


ca VG?) 50") 
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where the dressing functions S(p?) and V (p?) will be called the scalar and vector 
form factor respectively at below. The OPE of these two form factors renormalized 
in the MS scheme in Landau gauge was calculated to three loops in Ref. [16]. Up 
to operators of dimension three, one has 


Cms (fb, p? Cy ax (p, p 
Sap) = Sporlu pml) + SEP ys Coe P, y 
Clits p?) , - 
HERD gu), (2) 
and A g a r 
Vig?) = Vor + ERPI PRSE an (3) 


Here the purely perturbative parts Spr and Vpr were computed at three loops in 
Ref. [L7]. The Wilson coefficients Cms, Cj, 4», C5, Cm2 and C4» at three loops 
can be found in Ref. [16]. 

In principle if we can obtain the scalar and vector form factors by lattice QCD, 
then we can fit the lattice data to the functions in Eqs.(2]3) to extract out the 
quark mass and the chiral condensate. Since we need the inverse powers of p? to 
suppress the contributions from higher dimension operators, the lower limit of the 
fitting range in p? can not be too small. The Wilson coefficients are calculated 
by perturbation theory. This also requires p? can not be too small. On the other 
hand, if p? is too large then O(a?p?) and higher order lattice discretization effects 
in the data will be out of control. Thus one needs to find a fitting window in which 
a stable and reliable value for the chiral condensate can be obtained. 

Before the fittings we do not know if such a window exists or not given the 
lattice spacing in our data. Therefore we will vary our fitting range to test the 
reliability of our results. And we shall take into account the lattice discretization 
artifacts in our error analysis. 


3 Lattice setup 


We use the 2+1-flavor domain wall fermion configurations generated by the RBC- 
UKQCD collaborations [18]. The parameters of the ensembles used in this analysis 
are given in Tab. |1| Three light sea quark masses are used to check the sea quark 
mass dependence of our results. 

We use overlap fermions for the valence quark. The massless overlap opera- 
tor [20] is defined as 


Dov(p) = 14-458 Dp) (4) 
Here ¢ is the matrix sign function and D,y(p) is the usual Wilson fermion operator, 
except with a negative mass parameter —p = 1/2& — 4 in which ke < k < 0.25. In 
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Table 1: Parameters of the 2+1-flavor domain wall fermion configurations gener- 
ated by the RBC-UKQCD collaboration. The residual mass is from Ref. [18]. The 
lattice spacing was determined in Ref. [19]. 
l/a(GeV) label amj*"/am;^ volume Nu AMres 
1.75(4) | e005  0.005/0.04 24% x 64 92 0.003152(43) 
c01 0.01/0.04 — 24? x64 88 
c02 0.02/0.04 — 24? x 64 138 


Table 2: The overlap valence quark masses am, in lattice units used in this analysis. 
0.00620 0.00809 0.01020 0.01350 0.01720 0.02430 0.03650 0.04890 


our calculation we use & = 0.2 which corresponds to p = 1.5. The massive overlap 
Dirac operator is defined as 


Dm = pD«(p)*m(1— Pop) 
= p+” +(e- 2) q5e(1sDu(p)). (5) 
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To accommodate the SU(3) chiral transformation, it is usually convenient to use 
the chirally regulated field  — (1 — 5 Dov) in lieu of w in the interpolation field 
and operators. That is to say, our valence quark propagator is 


= Dy 
Ce), a= ; Dn = , (6) 


where De = m is chiral, i.e. (5s, De} = 0 [21]. 

The overlap valence quark masses in lattice units are given in Tab. The 
corresponding pion masses range from 220 to 600 MeV. Our quark propagators 
are calculated by using a point source on each configuration. The numbers of 
configurations used in this work are given in Tab. For three of the valence 
quark masses (0.01350, 0.02430, 0.04890) on ensemble c005, eight point sources 
on each configuration are used. For the same three quark masses on ensemble 
c02, eight point sources are used on half of the 138 configurations. The eight 
point sources are evenly distributed on the time slides and randomly distributed 
in 3-space from configuration to configuration to reduce autocorrelations. Part of 
these propagators were calculated and used in the computation of renormalization 
constants and in the study of diquarks [23]. We average the quark propagators 
from the eight sources on each configuration for these three valence quark masses. 
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Then together with the data from other configurations for other quark masses a 
Jackknife procedure (one configuration eliminated each time) is done to get the 
statistical uncertainties in our analysis below. Since ensemble c01 has the least 
statistics, the result from it will have the largest uncertainty. While c005 will have 
the smallest statistical uncertainty. 

Anti-periodic and periodic boundary conditions are used respectively in the 
time and spacial directions. Therefore the momentum modes are 


2nk, 2k, 2mka CL REDE (7) 
LT D XR — 


ap = ( 


where k,, = —6, —5,...,6. To reduce Lorentz noninvariant discretization effects, we 
use the momentum modes close to the diagonal line. This is achieved by doing a 
cut as was done in Ref. [22] 


[4] : 

p n n 

(pii < 0.32, where p”! = 1 p p= J p,. (8) 
p=1 H 


4 Analysis and discussions 


From Eq. (1) we have 


gr hisp] =, 


1 

1; Pris. ()] = V»). (9) 
In Fig. |1| we show the bare scalar and vector form factors (S(p?)/p? and V (p?)) 
in lattice units from our data ensemble c02 as functions of a?p?. The scalar form 
factor has a visible quark mass dependence as is shown in the graph on the left. On 
the contrary, the vector form factor in the graph on the right has no visible quark 
mass dependence even at quite low p? region. For example, at a?p? — 1.114 the 
vector form factors for am, — 0.00620 and 0.04890 agree with each other within the 
statistical uncertainties (0.720(4) versus 0.717(3)). This indicates the contribution 
from the m? term is quite small in Ba. Therefore we can also expect the m3 
term in Eq.(2) is negligible. Indeed at below we will see the quark mass dependence 
of the scalar form factor can be well described by a linear function. 

In our analysis below we take into account the reduced O(a?p?) discretization 
effects by adding a term proportional to a?p? in the fitting functions. However 
there are other artifacts of O(a?g?). In Ref. [I] the authors find that O(a?g?) 
effects are substantial in the vector form factor V, but modest in the ratio S/V. 
Since we have not computed the lattice artifacts of O(a?g?) and thus can not 
remove them from our form factors, we estimate their effects in our results by 
comparing the chiral condensates obtained from analyzing the ratio of the form 
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Figure 1: Left: The bare vector form factor S(p?) divided by a?p? from the quark 
propagators for various valence quark masses. Right: The bare vector form factor 
V (p?) for various valence quark masses. 


factors and from analyzing the scalar form factor alone. The difference in the 
results from the two analysis will be taken as a systematic uncertainty. 


4.1 Analysis of the ratio of scalar to vector form factor 


Since the ratio is expected to have much smaller O(a?g?) lattice artifacts than the 
scalar form factor, we trust more on the chiral condensate from the analysis of the 
ratio. The number from this analysis will be taken as our final result. 

The gluon condensate (A?) in Landau gauge was determined in, for exam- 
ple, Refs. [24] 25]. In the analysis of Ref. [11], a compatible value of (A?) was 
found but it seemed not yet stable against the order in perturbation theory. Since 
(A?) /p? <~ 0.6 GeV?/4 GeV? = 0.15 is small in the range of p? in our following 
analysis and the corresponding Wilson coefficient |C42| is also small (~ 0.3), we 
ignore the contribution from this condensate in Eq.(3] as a first step (note Vpr 
is of order 1). To obtain information about (A4?) from analyzing the vector form 
factor, we need more statistics and need to subtract the O(a?g?) artifacts. 

The quark mass dependence of the vector form factor is quite small as was seen 
in Fig. This indicates we can keep only the first term on the right hand side 
of Eq.(3] in analyzing our data. Thus from Eqs. (2]3). we have for the bare and 


renormalized form factors 
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Here the quark field renormalization constants Z, in the numerator and denomi- 
nator cancel each other. 
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Figure 2: Left: The ratio of form factors for various valence quark masses. Right: 
Examples of linear extrapolation to the quark massless limit at three typical mo- 
mentum values. 


Define a ratio " 
Tr[Ss(p)] S(p*) 


= = : 11 
Trips, DJ] - PV) E 
then in the chiral limit we have 
Crp) s 
lim R = A / , 12 
Im p)? Vpr (py) (12) 


In lattice units and taking into account O(a?p?) lattice artifacts in the quark 
propagator, we use the following function 
R = Cas (p, p?) 


. SI 2,2 
El a m (a2p2)2Vpp— wv) + Ba p. (13) 


to fit the ratio obtained from our lattice quark propagator. The dimensionless 
quark chiral condensate a? (iw) and B are two fit parameters. 

In the graph on the left of Fig. |2| we show the ratio R (in lattice units) as a 
function of a?p? from ensemble c02 for various valence quark masses. The graph 
on the right of Fig. |2| shows examples of the linear chiral extrapolation of R at 
three typical momentum values: a?p? = 1.024, 1.572 and 7.694. At all momentum 
values in our data for R we see a good linear dependence on amy. 

Then we fit the ratio R in the chiral limit to the function Eq. (13). Fig. [3]|shows 
an example of the fitting using a fitting range a?p? € [2.2, 5.3]. The fitting in the 
right graph does not include the Ba?p? term. Comparing it with the fitting in 
the left graph which does contain this term, we see that the Bap? term decreases 
x? /dof significantly. In evaluating the Wilson coefficients in the fitting function, we 
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Figure 3: Left: Fitting of the ratio R in the valence quark massless limit. Right: 
The same fit as in the left graph but without the Ba?p? term. The data points in 
purple (a?p? € [2.2,5.3]) are included in the fittings. 


use Ans — 332(17) MeV for three flavors in the MS scheme to compute the 


strong coupling constant œs. aMS(2 GeV) is calculated by using its perturbative 
running to 3-loops since the Wilson coefficients are only known to 3-loops. From 
this fitting, we get (uw)M5(2 GeV) = —(307(37)(7) MeV)? by using the lattice 
spacing 1/a — 1.75(4) GeV [19]. Here the first uncertainty is statistical and the 
second is from the uncertainty in the lattice spacing. 

To check the stability of the result against the fitting range, we vary the lower 
and upper limits of a?p?. In Tab. 3| we give the x?/dof of the fittings and the 
results of (aw) against these changes. As we see from the table, we can get a 
stable value for (uw). 

We then check the truncation error from the perturbative expansion of the 
Wilson coefficients. We repeat the fittings with Wilson coefficients and a, be- 
ing evaluated at 2-loops and 1-loop. The resulted numbers from data ensem- 
ble c02 are collected in Tab. Taking the difference between the center values 
with n = 2 and n = 3 as a systematic error, we finally get (i)M9(2 GeV) = 
(—307(37)(7)(8) MeV)? on ensemble c02. This is collected in Tab. 

Similarly in Tab. [5| and Tab. [6] we give the results from various fitting ranges 
on the other two ensembles c01 and c005 respectively. The truncation effects in 
the Wilson coefficients and a, are examined too. The quark condensates from all 
three ensembles are listed in Tab. 

We also tried to do fittings in a same momentum range a?p? € [2.0, 3.2] on all 
three ensembles. What we found are given in Tab. E 

Besides all the above, we repeat the fittings with (A?)MS(2 GeV) being fixed 
to 0.6 GeV? in Eq. (3). The resulted changes in ((iv))'? are 3 to 4 MeV, much 


Table 3: (ywW)™8(2 GeV) from fittings of R with different fitting ranges on ensemble 
c02. The first uncertainty is statistical and the second is from the uncertainty in 


the lattice spacing. 


Table 4: (uw)M9(2 GeV) from fittings of R on ensemble c02 with different trunca- 
tion order (n-loops) in evaluating the Wilson coefficients and as. The fitting range 


a^p € p'e/GeV? x?/dof ((ww))*/MeV 
2.2,5.5] [|6.7,16.8 0.91  —299(36)(7) 
2.2,5.3] [6.7162] | 0.71  —307(37)(7) 
22,51] [6.7,15.6] 0.73  —304(37)(7) 
2.24.9) [6.7,15.0] 0.76 ——305(38)(7) 
22,47] [|6.7,.144] 0.79 —309(39)(7) 
2.6,5.34 |8.0,16.2] 0.65  —299(41)(7) 
2.4,5.3]  [7.4,16.2] 0.61 | —300(39)(7) 
2.0,5.3) [6.1,16.2] 0.77  —310(35)(7) 
1.8,5.3] [5.5,16.2] 101  —318(32)(7) 


is a?p? € [2.2, 5.3]. 


Table 5: (uw)M9(2 GeV) from various fitting ranges on ensemble c01. The first 
uncertainty is statistical and the second is from the uncertainty in the lattice 
spacing. 


n x*/dof ((yw))/3/MeV 
1 0.73  —330(39)(8) 
2 072  —315(38)(7) 
3 Q7 -307(37)(7) 


ape  p^c/GeV? x?/dof ((ww))/3/MeV 
18,38) 55,116] 0.97 —299(42)(7) 
1.8,3.6] [5.5,11.0) 0.99  —300(43)(7) 
1.8,3.4] [5.5,10.4) 1.08  —293(45)(7) 
22,3.8 [6.7,11.6] 0.87 —304(53)(7) 
2.0,3.8] [6.1,1L6] 0.78 — —295(48)(7) 
1.6,3.8) [4.9,11.6] 131  —322(32)(7) 


Table 6: (uw)MS(2 GeV) from various fitting ranges on ensemble c005. The first 
uncertainty is statistical and the second is from the uncertainty in the lattice 
spacing. 


ap E p'e/GeV^ x^/dof ((ww))*/MeV 
14,3.] [43,101] 115  -302ü1)(7) 
14,3.] [4.3,9.5) 0.90  -306(11)(7) 
14,2.] [43,89] 0.998  —306(12)(7) 
L8,31]  |5.5,9. 0.81  —303(16)(7) 
16,3.1] [49,95] 0.73  —300(13)(7) 
1.2,3.1] ([3.75,9.5] 0.94  —310(10)(7) 


Table 7: (4) M9(2 GeV) on the three ensembles. The first uncertainty is statisti- 
cal. T'he second one is from the uncertainty of the lattice spacing. The third one is 
an estimation of the truncation effects in the evaluations of the Wilson coefficients 
and o;. 


ensemble  a?p?^€  p?€/GeV? x?/dof ((ww))'/3/MeV 
02 [253 [67163] 071  —307(37)(7)8) 
c0) —[L8,3.] [55,1100 0.97  —299(42)(7)(11) 
c005  [L431] [43,95] 090  -306(11)(7)(13) 


Table 8: (uw)M5(2 GeV) from a same fitting window a?p? € [2.0, 3.2] on all three 
ensembles. The three uncertainties are as explained in Tab. 
ensemble x?/dof ((ww))'?/MeV 
c02 0.78  —330(39)(8)12) 
c01 0.79 —281(57)(6)(10) 
c005 0.85  —303(18)(7)(12) 
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Figure 4: The light sea quark mass dependence of (yw)M5(2 GeV). The crosses 
are from constant fits. 


smaller than the statistical or other uncertainties. This means it is safe to ignore 
the contribution from (A?) with our current statistics. 

The light sea quark mass dependence is shown in Fig. 4| where we plot together 
the results from all three ensembles. The three red points are those in Tab. [7] from 
fittings with different p? range on each ensemble. The blue ones are those in Tab. 
from fittings in a same momentum range on all three ensembles. In this graph, we 
have quadratically combined together the three uncertainties in Tab. [7| and Tab. 
respectively. Since we do not see an apparent sea quark mass dependence with our 
relatively large uncertainties, we do a constant fit to finally obtain 


(bw) 5 (2 GeV) = (—305(15) MeV)? (different fitting ranges), (14) 


and 
(b)M$(2 GeV) = (—306(19) MeV)? (same fitting range). (15) 


These two numbers are in good agreement with each other. 


4.2 Analysis of the Scalar form factor 


There may be non-negligible O (a?g?) lattice artifact in our scalar and vector form 
factors as were seen in Ref. [IT] with Wilson twisted mass fermions. At large 
p’, difference was seen in the O(a?g?)-corrected and un-corrected vector form fac- 
tor [11]. Unfortunately, We have not calculated this artifact yet and therefore 
could not do this correction to our data. To estimate its effects, we analyze the 
scalar form factor in the chiral limit to obtain the chiral condensate and compare 
the result with the one from Sec. 
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Figure 5: The quark field renormalization constant for ensemble c005. The black 
vertical line indicates the position of u = 2 GeV. 


From Eq.(2) we see in the chiral limit the scalar form factor is related to the 
chiral condensate by 


Coy (Hs p?) 
Sale") = 09) Q0. (16) 
With a quark field renormalization constant Yr = zy ^a and taking into account 
O(a*p*) effects, we have 


Se(a^p?)  Sa(a?p?) _ Cy (m, p^) 
ap? ZQa2p? Z,(o?p^y 


rS, (a, p)] = 


12 a (pp) (u) + Bap’. (17) 


Here we have put everything in lattice units. Thus the quark field renormalization 
constant Z, is needed in the analysis of the scalar form factor. 


4.2.1 Quark field renormalization 


Our Z, is first calculated in the RI-MOM scheme and then converted to the 
MS scheme. The detailed calculation in the RI-MOM scheme for our work can 
be found in Ref. [22]. We first use the axial vector Ward Identity to obtain Z/, 
which equals to Z4 in the RI-MOM scheme. Then from it Z, in the RI-MOM 
scheme is computed at several valence quark masses. The results of Z, show little 
quark mass dependence (see Fig.3 in Ref. [22]). We now do a linear extrapolation 
of Z, in the quark mass to the chiral limit. The results in this limit are shown by 
the black diamonds in Fig.|5| Then the conversion ratio calculated by perturbation 
theory [I7] to 3-loops is used to get Z, in the MS scheme, which is shown by the 
red crosses in Pi B After running ZM’ from an initial scale p? to u? = (2 GeV)? 
by using its anomalous dimension to 3-loops, we obtain the blue pluses in Fig. 
The deviation of the blue pluses from a constant at large initial scales is attributed 
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Table 9: ZN (2 GeV) on the three ensembles. The first error is statistical and the 
second one is a 1% systematic error. 
ensemble c02 c01 c005 
ZMS(2 GeV) 1.202(2)(12) 1.209(3)(12) 1.197(2)(12) 
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Figure 6: Left: Examples of linear extrapolation of the scalar form factor to the 
quark massless limit at three typical momentum values. Right: Fitting of the 
scalar form factor in the chiral limit to Eq.( on ensemble c02. 


to O(a?y?) lattice artifacts. Thus a linear extrapolation in a?p? to a?p? = 0 is done 
to get p. GeV) (illustrated by the blue line using data points at a?p? > 5). 

The results of Z, in the MS scheme are given in Tab. [9| for the three ensembles. 
Similarly to what have been done to Zs for the scalar density in Ref. (see 
its Tab. V), we find a 1% systematic uncertainty for Z, from the uncertainty in 
the lattice spacing, the uncertainty in APO the truncation error of the pertur- 
bative conversion ratio and the variation of the fitting range in the a?p? linear 
extrapolation. This systematic uncertainty is given in Tab. [9] 


4.2.2 Fitting results 


We have shown the scalar form factor (divided by p?) in the left graph of Fig. 
In Fig. [6] we show the linear chiral extrapolation of the scalar form factor and the 
fitting of the chiral limit results to the function Eq.(17]. Again, we find in the fit 
the Bap? term decreases X?/dof significantly. As was done in the analysis of the 
ratio of form factors in Sec. we check the stability of the results of (ii) against 
the fitting range in p”, and against the order of truncation in the evaluations of the 
Wilson coefficients and œs. Since the uncertainty of Z; is quite small compared with 
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Table 10: (4w)M9(2 GeV) on the three ensembles from the analysis of the scalar 
form factor. The first uncertainty is statistical. The second one is from the un- 
certainty of the lattice spacing. The third one is an estimation of the truncation 
effects in the evaluations of the Wilson coefficients and as. 
ensemble  a?p^€  p?€/GeV? x?/dof ((ww))'/8/MeV 

c02 [2.2, 5.3] — [6.7, 16.2] 0.706 | —272(32) 

c01 [1.8,3.8]  [5.5, 11.6] 1.00  —278(35) 

c005 — [L4,3.] [43,9.5] 103  —288(10) 


Table 11: (uw)M9(2 GeV) from fittings of the scalar form factor in a same fitting 
window a?p? € [2.0,3.2] on all three ensembles. The three uncertainties are as 
explained in Tab. 


ensemble x?/dof ((ww))/3/MeV 
c02 1.00 —299(35)(7)(9) 
c01 0.83 —262(47)(6)(9) 
c005 0.86 —280(16)(6)(9) 


other sources of uncertainties, we have ignored its propagation to the uncertainty 
of the quark chiral condensate. 

The results of (ww) from the three ensembles are given in Tabs. Tab 
is for fittings with different p? window on each ensemble and Tab |11|for fittings 
with a same p? window on all three ensembles. They are in agreement within 
errors. 

The light sea quark mass dependence of (Ww) is again small compared with 
our uncertainties. Thus we do a constant fit similar to Fig. 4| in Sec. What 
we obtain are 


(bw) M9(2 GeV) = (—284(13) MeV)? (different fitting ranges), (18) 


and 


Qu) (2 GeV) = (—282(16) MeV)? (same fitting range). (19) 


5 Summary 
We determine the quark chiral condensate by fitting lattice data of the overlap 


quark propagator to its operator product expansion in the MS scheme in Landau 
gauge. We perform two analyses. One uses the ratio of scalar to vector form factor 
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of the propagator, which is supposed to have modest O(a?g?) lattice artifacts. The 
other one uses the scalar form factor. We use the result from the second analysis 
to estimate the uncertainty from the O(a?g?) artifacts. The fitting range of the 
momentum in our analysis is varied to check the stability of the results. The 
truncation error in evaluating the Wilson coefficients and a, is also examined. 
Three ensembles of 2+1-flavor domain wall fermion configurations are used to 
check the light sea quark mass dependence. 

We take the number in Eq.(14) as our final result. The difference between the 
center values in Eq.(14) and ME is taken as a systematic uncertainty due to 
the O(a?g?) effects in our data. That is to say, our final result is 


(bw) MS(2 GeV) = (—305(15)(21) MeV)?. (20) 


Here the first error contains uncertainties from statistics, the lattice spacing and 
truncations in perturbative calculations of the Wilson coefficients and a,. 

Our result Eq. (20), with a relatively large error bar, agrees with the FLAG-3 
average 41/3 = 274(3) MeV for Nr = 2 + 1 flavor lattice calculations. To 
improve our work, the O(a?g?) effects should be calculated and removed from the 
lattice data of the quark propagator. With more statistics the effects of the (A?) 
term can be checked carefully. Furthermore, calculations at more lattice spacings 
should be done to enable a continuum extrapolation. 
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